unit NumerousMethods;

{
  PasCalc
  Version 1.2 2011-05-26

  Copyright (C) 2011 KEEPER <yar.keeper@gmail.com>

  This file is part of the PasCalc Calculator

  This unit is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.

  See the file LICENSE.TXT, included in this distribution,
  for details about the copyright.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
}

{
  ×èñëåííûå ìåòîäû
  (c) Keeper
  (r) Âàñèëüêîâà Í. Í, Âàñèëüêîâ Þ. Â.
}

interface

uses
  Math, MatrixOperations, Analyser, Graphics, SysUtils, ComplexOps, Classes;

type

  //Òàáëè÷íî çàäàííàÿ ôóíêöèÿ âèäà
  //|x0|x1|..|xn|
  //|y0|y1|..|yn|
  TFunction = array [0..1] of TVector;  //0 - args; 1 - values

  function CopyFunction(const Source: TFunction): TFunction;

type

  TPoint = record
    x, y: integer;
  end;

  TFloatPoint = record
    x, y: extended;
  end;

  (*
    Òàê êàê ïðîäîëæèòåëüíîñòü ïåðåõîäíîãî ïðîöåññà äî íà÷àëà
    âû÷èñëåíèé íåèçâåñòíà, à äëèíà òðóáû äåòåðìèíèðîâàííà, òî
    îïèøåì äèíàìèêó êàê äèíàìè÷åñêèé ìàññèâ ñïèñêîâ
  *)

  //Ñòðóêòóðà, îïèñûâàþùàÿ òî÷êó äèíàìèêè
  TDynPoint = record
    Time, Value: extended
  end;

  PDynPoint = ^TDynPoint;

  //Äèíàìèêà ïðîöåññà íà ó÷àñòêå
  TProcessRegion = class(TList)
  private
    function Get(Index: Integer): TDynPoint;
    procedure Put(Index: Integer; const Value: TDynPoint);
  public
    Distance: extended;
    function Add(Item: TDynPoint): Integer;
    procedure Delete(Index: Integer);
    property Items[Index: Integer]: TDynPoint read Get write Put; default;
  end;

  { Îáîáù¸ííûé èíòåðïîëÿòîð. Ñîäåðæèò àáñòðàêòíûå ìåòîäû }
  TInterpolator = class
  private
    fSourcePoints: TFunction; //Èñõîäíàÿ òàáëèöà çíà÷åíèé ôóíêöèè
    fInterpolant: TVector; //Âåêòîð ñ ðåçóëüòàòîì èíòåðïîëÿöèè
    fUnknownPoints: TVector; //Âåêòîð ñ íåèçâåñòíûìè òî÷êàìè
    procedure SetSourcePoints(const AValue: TFunction);
    procedure SetUnknownPoints(const AValue: TVector);
    //Ïîäðàçóìåâàåòñÿ, ÷òî áîëüøàÿ ÷àñòü íåîáõîäèìûõ âû÷èñëåíèé
    //âûïîëíÿåòñÿ ïîñëå çàäàíèÿ òàáëèöû çíà÷åíèé ôóíêöèè. Äàííûé ìåòîä
    //ÂÑÅÃÄÀ äîëæåí áûòü ïåðåîïðåäåë¸í!
    procedure AfterSetSourcePoints; virtual; abstract;
  public
    { Êîíñòðóêòîðó ïåðåäà¸òñÿ òàáëèöà çíà÷åíèé ôóíêöèè SourceTable }
    constructor Create(SourceTable: TFunction);
    { Òàáëèöà ñî çíà÷åíèÿìè ôóíêöèè }
    property SourcePoints: TFunction read fSourcePoints write SetSourcePoints;
    { Ìàññèâ ñ íåèçâåñòíûìè òî÷êàìè }
    property UnknownPoints: TVector read fUnknownPoints write SetUnknownPoints;
    { Ìàññèâ ñ ðåçóëüòàòîì èíòåðïîëÿöèè }
    property InterpResult: TVector read fInterpolant;
    { Çàïîëíÿåò âåêòîð fInterpolant ðàññ÷èòàííûìè ïî fUnknownPoints çíà÷åíèÿìè }
    procedure FindUnknownPoints;
    { Ìåòîä ïîèñêà çíà÷åíèÿ â êîíêðåòíîé òî÷êå Õ. Óíèêàëåí äëÿ êîíêðåòíîãî
      èíòåðïîëÿòîðà. Âñåãäà äîëæåí áûòü ïåðåîïðåäåëåí! }
    function GetPoint(const X: extended): Extended; virtual; abstract;
  end;

  { Èíòåðïîëÿòîð ìåòîäîì Ëàãðàíæà }
  TLagrangeInterpolator = class(TInterpolator)
  private
    fQDenominators: TVector;
    fQVector: TVector;
    procedure FormQDenominatorVector;
    function GetQDenominator(const idx: integer): extended;
    function GetQNumenator(const idx: integer; const X: extended): extended;
    procedure AfterSetSourcePoints; override;
    function GetQPoint(i: integer): extended;
  public
    property GetQ[i: integer]: extended read GetQPoint;
    function GetPoint(const X: extended): Extended; override;
    function QCount: integer;
  end;

  TNewtonInterpolator = class(TInterpolator)
  private
    fYDeltas: TVector;
    fHValue: extended;
    fXDeltas: TVector;
    Matrix: TMatrix;
    procedure CalcH;
    procedure CalcXDeltas(const X: extended);
    procedure CreateYDeltas;
    procedure AfterSetSourcePoints; override;
  public
    property YDeltas: TMatrix read Matrix;
    property H: extended read fHValue;
    function GetPoint(const X: extended): Extended; override;
  end;

  TApproximator = class(TInterpolator)
  private
    fDegree: integer;
    fDegreesVector: TVector;
    EquationSystem: TMatrix;
    fSolutionVector: TVector;
    fXYSum: TVector;
    procedure FormXYSumVector;
    procedure FormEquationSystem;
    procedure AfterSetSourcePoints; override;
    procedure GetDegreesVector;
    procedure SetDegree(const Value: integer);
    function GetPolynomCoefficient(i: integer): extended;
  public
    constructor Create(SourceTable: TFunction; PolynomDegree: integer);
    property Matr: TMatrix read EquationSystem;
    property Degrees: TVector read fDegreesVector;
    property XYsum: TVector read fXYSum;
    property PolynomCoefficient[i: integer]: extended read GetPolynomCoefficient;
    function CoefficientsCount: integer;
    property Degree: integer read fDegree write SetDegree;
    function GetPoint(const X: extended): Extended; override;
  end;

  { Èíòåãðèðîâàíèå ïðè ïîìîùè êâàäðàòóðíûõ ôîðìóë - îáîáù¸ííûé êëàññ }

  { TQuadratureIntegrator }

  TQuadratureIntegrator = class
  private
    fDegree: integer;
    fIsRad: boolean;
    fLeft, fRight: extended;
    fArguments, fValues: TVector;
    fFunction: string;
    fCalculator: TCalculator;
    fNewtonKotesCoefs: TVector;
    fMiddleSums: TVector;
    fDivide: integer;
    fWi: TVector;
    fIntegrationResult: extended;
    procedure GetNewtonKotesCoefs;
    procedure CreateArgumentsVector;
    procedure CalcMiddleSums;
    procedure SetDegree(const AValue: integer);
    procedure SetIsRad(const AValue: boolean);
    procedure CalcWi;
    procedure SetFunction(const Value: string);
    procedure SetLeft(const Value: extended);
    procedure SetRightt(const Value: extended);
  public
    { Ïîëó÷àåò ïðåäñòàâëåíèå óãëîâ è ñòåïåíü ïîëèíîìà }
    constructor Create(const IsRad: boolean; const Degree: integer);
    { Èíòåãðèðîâàòü ôóíêöèþ Func }
    function Integrate: extended;
    { Ñòåïåíü ïîëèíîìà (1..5). Óâåëè÷åíèå ïîâûøàåò êà÷åñòâî è ñíèæàåò ñêîðîñòü }
    property PolynomDegree: integer read fDegree write SetDegree;
    { Ïðåäñòàâëåíèå óãëîâ }
    property IsRadians: boolean read fIsRad write SetIsRad;
    { Èíòåãðèðóåìàÿ ôóíêöèÿ }
    property Func: string read fFunction write SetFunction;
    { Ëåâûé ïðåäåë èíòåãðèðîâàíèÿ }
    property LeftLimit: extended read fLeft write SetLeft;
    { Ïðàâûé ïðåäåë èíòåãðèðîâàíèÿ }
    property RightLimit: extended read fRight write SetRightt;
  end;

  TSolverType = (stXY, stY);

  TRungeKuttSolver = class
  private
    fStartCond: TFloatPoint;
    fStep: extended;
    fFunction: string;
    fFinishParam: extended;
    fSolveResult: TFunction;
    fSolveToResult: TProcessRegion;
    fSolverType: TSolverType;
    procedure SetStartCond(const Value: TFloatPoint);
    procedure SetStep(const Value: extended);
    procedure SetFunction(const Value: string);
    procedure SetFinishParam(const Value: extended);
    procedure CalculatePointsNormal;
    procedure SeSolverType(const Value: TSolverType);
  public
    destructor Destroy; override;
    property StartCond: TFloatPoint read fStartCond write SetStartCond;
    property Step: extended read fStep write SetStep;
    property RFunction: string read fFunction write SetFunction;
    property FinishParam: extended read fFinishParam write SetFinishParam;
    { Åñëè ïðàâàÿ ÷àñòü ôóíêöèÿ ïåðåìåííîé è ðåøåíèÿ }
    procedure Solve;
    procedure SolveTo(const Limit: extended);
//    function GetPoint: extended;
    { Åñëè ïðàâàÿ ÷àñòü ôóíêöèÿ ðåøåíèÿ }
    property SolveResult: TFunction read fSolveResult;
    property SolveToResult: TProcessRegion read fSolveToResult;
    property SolverType: TSolverType read fSolverType write SeSolverType;
  end;

implementation

{ TQuadratureIntegrator }

procedure TQuadratureIntegrator.GetNewtonKotesCoefs;
begin
  {
    Îãðàíè÷èâàåì
  }
  if fDegree > 5 then fDegree := 5;
  if fDegree < 1 then fDegree := 1;
  SetLength(fNewtonKotesCoefs, fDegree + 1);
  case fDegree of
    1:
    begin
      fNewtonKotesCoefs[0] := 0.5;
      fNewtonKotesCoefs[1] := 0.5;
    end;
    2:
    begin
      fNewtonKotesCoefs[0] := 1/6;
      fNewtonKotesCoefs[1] := 2/3;
      fNewtonKotesCoefs[2] := 1/6;
    end;
    3:
    begin
      fNewtonKotesCoefs[0] := 1/8;
      fNewtonKotesCoefs[1] := 3/8;
      fNewtonKotesCoefs[2] := 3/8;
      fNewtonKotesCoefs[3] := 1/8;
    end;
    4:
    begin
      fNewtonKotesCoefs[0] := 7/90;
      fNewtonKotesCoefs[1] := 16/45;
      fNewtonKotesCoefs[2] := 2/15;
      fNewtonKotesCoefs[3] := 16/45;
      fNewtonKotesCoefs[4] := 7/90;
    end;
    5:
    begin
      fNewtonKotesCoefs[0] := 19/288;
      fNewtonKotesCoefs[1] := 25/96;
      fNewtonKotesCoefs[2] := 25/144;
      fNewtonKotesCoefs[3] := 25/144;
      fNewtonKotesCoefs[4] := 25/96;
      fNewtonKotesCoefs[5] := 19/288;
    end;
  end;
end;

procedure TQuadratureIntegrator.CreateArgumentsVector;
var
  PointsCount, i: integer;
  Step: extended;
  Value: extended;
begin
  PointsCount := (fDivide * (fDegree)) + 1;
  SetLength(fArguments, PointsCount);
  Step := (fRight - fLeft)/(PointsCount - 1);
  Value := fLeft;
  for i := 0 to High(fArguments) do
  begin
    fArguments[i] := Value;
    Value := Value + Step;
  end;
end;

procedure TQuadratureIntegrator.CalcMiddleSums;
var
  i, j: integer;
  Cursor: integer;
  test: extended;
begin
  fCalculator := TCalculator.Create(fFunction);
  fCalculator.IsRad := fIsRad;
  fValues := fCalculator.GetValuesFromArguments(fArguments);
  fCalculator.Free;
  SetLength(fMiddleSums, fDivide);
  Cursor := 0;
  for i := 0 to High(fMiddleSums) do
  begin
    fMiddleSums[i] := 0;
    for j := 0 to High(fWi) do
    begin
      Test := fValues[cursor] * fWi[j];
      fMiddleSums[i] := fMiddleSums[i] + test;
      inc(cursor);
    end;
    dec(cursor);
  end;
end;

procedure TQuadratureIntegrator.SetDegree(const AValue: integer);
begin
  if fDegree = AValue then exit;
  fDegree := AValue;
end;

procedure TQuadratureIntegrator.SetIsRad(const AValue: boolean);
begin
  if fIsRad = AValue then exit;
  fIsRad := AValue;
end;

constructor TQuadratureIntegrator.Create(const IsRad: boolean;
  const Degree: integer);
const
  ACCURACY = 200; //×èñëî ïîäûíòåðâàëîâ. ×åì áîëüøå òåì âûøå òî÷íîñòü è ìåäëåííåé ðàáîòà
begin
  fDivide := ACCURACY;
  fIsRad := IsRad;
  fDegree := Degree;
  GetNewtonKotesCoefs;
end;

function TQuadratureIntegrator.Integrate: extended;
var
  i: integer;
begin
  CreateArgumentsVector;
  CalcWi;
  CalcMiddleSums;
  result := 0;
  for i := 0 to high(fMiddleSums) do
    result := result + fMiddleSums[i];
  fIntegrationResult := result;
end;

procedure TQuadratureIntegrator.CalcWi;
var
  bigstep: extended;
  i: integer;
begin
  bigstep := (fRight-fLeft)/fDivide;
  SetLength(fWi, Length(fNewtonKotesCoefs));
  for i := 0 to High(fNewtonKotesCoefs) do
    fWi[i] := bigstep * fNewtonKotesCoefs[i];
end;

procedure TQuadratureIntegrator.SetFunction(const Value: string);
begin
  fFunction := Value;
end;

procedure TQuadratureIntegrator.SetLeft(const Value: extended);
begin
  fLeft := Value;
end;

procedure TQuadratureIntegrator.SetRightt(const Value: extended);
begin
  fRight := Value;
end;

{ TLagrangeInterpolator }

procedure TLagrangeInterpolator.AfterSetSourcePoints;
begin
  FormQDenominatorVector;
end;

procedure TLagrangeInterpolator.FormQDenominatorVector;
var
  i: integer;
begin
  SetLength(fQDenominators, Length(fSourcePoints[0]));
  for i := 0 to high(fQDenominators) do
    fQDenominators[i] := GetQDenominator(i);
end;

constructor TInterpolator.Create(SourceTable: TFunction);
begin
  SetSourcePoints(SourceTable);
end;

procedure TInterpolator.FindUnknownPoints;
var
  i: integer;
begin
  SetLength(fInterpolant, Length(fUnknownPoints));
  for I := 0 to High(fUnknownPoints) do
    fInterpolant[i] := GetPoint(fUnknownPoints[i]);
end;

procedure TInterpolator.SetSourcePoints(const AValue: TFunction);
begin
  fSourcePoints := AValue;
  AfterSetSourcePoints;
end;

procedure TInterpolator.SetUnknownPoints(const AValue: TVector);
begin
  fUnknownPoints := AValue;
end;

function TLagrangeInterpolator.GetQDenominator(const idx: integer): extended;
var
  i: integer;
begin
  result := 1;
  for i := 0 to idx - 1 do
    result := result * (fSourcePoints[0, idx] - fSourcePoints[0, i]);
  for i := idx + 1 to high(fSourcePoints[0]) do
    result := result * (fSourcePoints[0, idx] - fSourcePoints[0, i]);
end;

function TLagrangeInterpolator.GetQNumenator(const idx: integer;
  const X: extended): extended;
var
  i: integer;
begin
  result := 1;
  for i := 0 to idx - 1 do
    result := result * (X - fSourcePoints[0, i]);
  for i := idx + 1 to high(fSourcePoints[0]) do
    result := result * (X - fSourcePoints[0, i]);
end;

function TLagrangeInterpolator.GetPoint(const X: extended): Extended;
var
  i: integer;
begin
  result := 0;
  SetLength(fQVector, Length(fQDenominators));
  for i := 0 to high(fSourcePoints[0]) do
  begin
    fQVector[i] := (GetQNumenator(i, X))/(fQDenominators[i]);
    result := result + (fSourcePoints[1, i] * (fQVector[i]));
  end;
end;

function TLagrangeInterpolator.QCount: integer;
begin
  result := Length(fQVector);
end;

function TLagrangeInterpolator.GetQPoint(i: integer): extended;
begin
  result := fQVector[i];
end;

{ TNewtonInterpolator }

procedure TNewtonInterpolator.AfterSetSourcePoints;
begin
  CreateYDeltas;
  CalcH;
end;

procedure TNewtonInterpolator.CalcH;
begin
  fHValue := fSourcePoints[0, 1] - fSourcePoints[0, 0];
end;

procedure TNewtonInterpolator.CalcXDeltas(const X: extended);
var
  i: integer;
  Val: extended;
begin
  SetLength(fXDeltas, High(fSourcePoints[0]));
  val := 1;
  for I := 0 to High(fXDeltas) do
  begin
    val := val * (X - fSourcePoints[0, i]);
    fXdeltas[i] := val;
  end;
end;

procedure TNewtonInterpolator.CreateYDeltas;
var
  i, j: integer;
begin
  SetLength(Matrix, Length(fSourcePoints[0]), Length(fSourcePoints[0]));
  for I := 0 to high(Matrix) do
    Matrix[i, 0] := fSourcePoints[1, i];
  for i := 1 to High(Matrix) do
    for j := 0 to High(Matrix) - i do
      Matrix[j, i] := Matrix[j+1, i-1] - Matrix[j, i-1];
  SetLength(fYDeltas, High(Matrix));
  for I := 1 to Length(fYDeltas) do
    fYDeltas[i - 1] := Matrix[0, i];
end;

function TNewtonInterpolator.GetPoint(const X: extended): Extended;

function factorial(const n: integer): extended;
var
  i: integer;
begin
  result := 1;
  for I := 2 to n do
    result := result * i;
end;

var
  I: Integer;
  Numenator, Denominator: extended;
begin
  CalcXDeltas(X);
  result := fSourcePoints[1, 0];
  for I := 1 to High(fSourcePoints[0]) do
  begin
    Numenator := fYDeltas[i-1]*fXDeltas[i-1];
    Denominator := factorial(i)*power(fHValue, i);
    result := result + Numenator/Denominator;
  end;
end;

{ TApproximator }

procedure TApproximator.AfterSetSourcePoints;
begin
  FormEquationSystem;
  fSolutionVector := TVector(GausseMethod(EquationSystem));
end;

function TApproximator.CoefficientsCount: integer;
begin
  result := Length(fSolutionVector);
end;

constructor TApproximator.Create(SourceTable: TFunction; PolynomDegree: integer);
begin
  fDegree := PolynomDegree;
  inherited Create(SourceTable);
end;

procedure TApproximator.FormEquationSystem;
var
  i, j: integer;
begin
  SetLength(EquationSystem, Degree + 1, Degree + 2);
  GetDegreesVector;
  for i := 0 to high(EquationSystem) do
    for j := 0 to high(EquationSystem[i]) - 1 do
      EquationSystem[i, j] := fDegreesVector[j + i];
  FormXYSumVector;
  for i := 0 to high(EquationSystem) do
    EquationSystem[i, high(EquationSystem[0])] := fXYSum[i];
end;

procedure TApproximator.FormXYSumVector;
var
  i, j: integer;
  sum: extended;
begin
  SetLength(fXYSum, Length(EquationSystem));
  for i := 0 to high(fXYSum) do
  begin
    sum := 0;
    for j := 0 to high(fSourcePoints[0]) do
      sum := sum + (fSourcePoints[1, j] * Power(fSourcePoints[0, j], i));
    fXYSum[i] := sum;
  end;
end;

procedure TApproximator.GetDegreesVector;
var
  i, j: integer;
  sum: extended;
begin
  SetLength(fDegreesVector, 2*Degree + 1);
  for i := 0 to high(fDegreesVector) do
  begin
    sum := 0;
    for j := 0 to high(FSourcePoints[0]) do
      sum := sum + power(fSourcePoints[0, j], i);
      fDegreesVector[i] := sum;
  end;
end;

function TApproximator.GetPoint(const X: extended): Extended;
var
  i: integer;
begin
  result := 0;
  for i := 0 to high(fSolutionVector) do
    result := result + fSolutionVector[i] * power(X, i);
end;

function TApproximator.GetPolynomCoefficient(i: integer): extended;
begin
  result := fSolutionVector[i];
end;

procedure TApproximator.SetDegree(const Value: integer);
begin
  fDegree := Value;
end;

{ TRungeKuttSolver }

procedure TRungeKuttSolver.CalculatePointsNormal;
var
  i: integer;
  k1, k2, k3, k4: extended;
  LexAnalyser: TLexAnalyser;
  SyntaxAnalyser: TSyntaxAnalyser;
  LexSeq: TLexemeArray;
  Variables: TComplexArr;
  err: TErr;
begin
  if fSolverType = stXY then
    SetLength(Variables, 2)
  else
    SetLength(Variables, 1);
  //Çàïóñêàåì ëåêñè÷åñêèé àíàëèçàòîð
  LexAnalyser := TLexAnalyser.Create(fFunction);
  try
    LexAnalyser.IsExpMode := false;
    LexAnalyser.Analyse(LexSeq); //Ðåçóëüòàò àíàëèçà â ïåðåìåííîé LexSeq
  finally
    LexAnalyser.Free;
  end;
  //Çàïóñêàåì ñèíòàêñè÷åñêèé àíàëèçàòîð
  SyntaxAnalyser := TSyntaxAnalyser.Create(LexSeq);
  try
    for i := 1 to high(fSolveResult[0]) do
    begin
      //Åñëè â ïðàâîé ÷àñòè 2 ïåðåìåííûå òî ìàññèâ èç 2 ýëåìåíòîâ
      if fSolverType = stXY then
      begin
        Variables[0].re := fSolveResult[0, i-1]; Variables[0].im := 0;
        Variables[1].re := fSolveResult[1, i-1]; Variables[1].im := 0;
      end else
      begin
        Variables[0].re := fSolveResult[1, i-1]; Variables[0].im := 0;
      end;
      //Ïåðåìåííàÿ óñòàíàâëèâàåòñÿ ÎÄÈÍ ðàç, äàëåå îíà ìåíÿåòñÿ àâòîìàòè÷åñêè,
      //òàê êàê äèíàìè÷åñêèé ìàññèâ ýòî óêàçàòåëü
      SyntaxAnalyser.Variable := Variables;
      k1 := fStep * SyntaxAnalyser.calc(err).re;
      if fSolverType = stXY then
      begin
        Variables[0].re := fSolveResult[0, i-1] + 0.5*step; Variables[0].im := 0;
        Variables[1].re := fSolveResult[1, i-1] + 0.5*k1; Variables[1].im := 0;
      end else
      begin
        Variables[0].re := fSolveResult[1, i-1] + 0.5*k1; Variables[0].im := 0;
      end;
      k2 := fStep* SyntaxAnalyser.calc(err).re;
      if fSolverType = stXY then
      begin
        Variables[0].re := fSolveResult[0, i-1] + 0.5*step; Variables[0].im := 0;
        Variables[1].re := fSolveResult[1, i-1] + 0.5*k2; Variables[1].im := 0;
      end else
      begin
        Variables[0].re := fSolveResult[1, i-1] + 0.5*k1; Variables[0].im := 0;
      end;
      k3 := fStep* SyntaxAnalyser.calc(err).re;
      if fSolverType = stXY then
      begin
        Variables[0].re := fSolveResult[0, i-1] + step; Variables[0].im := 0;
        Variables[1].re := fSolveResult[1, i-1] + k3; Variables[1].im := 0;
      end else
      begin
        Variables[0].re := fSolveResult[1, i-1] + k3; Variables[0].im := 0;
      end;
      k4 := fStep* SyntaxAnalyser.calc(err).re;
      fSolveResult[1, i] := fSolveResult[1, i-1] + ((k1+2*k2+2*k3+k4)/6);
    end;
  finally
    SyntaxAnalyser.Free;
  end;
end;


destructor TRungeKuttSolver.Destroy;
begin
end;

procedure TRungeKuttSolver.SeSolverType(const Value: TSolverType);
begin
  fSolverType := Value;
end;

procedure TRungeKuttSolver.SetFinishParam(const Value: extended);
begin
  fFinishParam := Value;
end;

procedure TRungeKuttSolver.SetFunction(const Value: string);
begin
  fFunction := Value;
end;

procedure TRungeKuttSolver.SetStartCond(const Value: TFloatPoint);
begin
  fStartCond.x := Value.x;
  fStartCond.y := Value.y;
end;

procedure TRungeKuttSolver.SetStep(const Value: extended);
begin
  fStep := Value;
end;

procedure TRungeKuttSolver.Solve;
var
  PointsCount: integer;
  i: integer;
begin
  if fFunction = '' then exit;
  PointsCount := Round((fFinishParam - fStartCond.x)/fStep)+1;
  SetLength(fSolveResult[0], PointsCount);
  SetLength(fSolveResult[1], PointsCount);
  fSolveResult[0, 0] := fStartCond.x;
  fSolveResult[1, 0] := fStartCond.y;
  for i := 1 to High(fSolveResult[0]) do
    fSolveResult[0, i] := fSolveResult[0, i-1] + fStep;
  CalculatePointsNormal;
end;


procedure TRungeKuttSolver.SolveTo(const Limit: extended);
var
  i: integer;
  Variables: TComplexArr;
  LexAnalyser: TLexAnalyser;
  LexSeq: TLexemeArray;
  SyntaxAnalyser: TSyntaxAnalyser;
  k1, k2, k3, k4: extended;
  err: Terr;
  Point: TDynPoint;
begin
  fSolveToResult := TProcessRegion.Create;
  Point.Time := fStartCond.x;
  Point.Value := fStartCond.y;
  fSolveToResult.Add(Point);

  if fSolverType = stXY then
    SetLength(Variables, 2)
  else
    SetLength(Variables, 1);
  //Çàïóñêàåì ëåêñè÷åñêèé àíàëèçàòîð
  LexAnalyser := TLexAnalyser.Create(fFunction);
  try
    LexAnalyser.IsExpMode := false;
    LexAnalyser.Analyse(LexSeq); //Ðåçóëüòàò àíàëèçà â ïåðåìåííîé LexSeq
  finally
    LexAnalyser.Free;
  end;
  //Çàïóñêàåì ñèíòàêñè÷åñêèé àíàëèçàòîð
  SyntaxAnalyser := TSyntaxAnalyser.Create(LexSeq);
  try
    i := 1;
    repeat
      //Åñëè â ïðàâîé ÷àñòè 2 ïåðåìåííûå òî ìàññèâ èç 2 ýëåìåíòîâ
      if fSolverType = stXY then
      begin
        Variables[0].re := fSolveToResult.Items[i-1].Time; Variables[0].im := 0;
        Variables[1].re := fSolveToResult.Items[i-1].Value; Variables[1].im := 0;
      end else
      begin
        Variables[0].re := fSolveToResult.Items[i-1].Value;; Variables[0].im := 0;
      end;
      //Ïåðåìåííàÿ óñòàíàâëèâàåòñÿ ÎÄÈÍ ðàç, äàëåå îíà ìåíÿåòñÿ àâòîìàòè÷åñêè,
      //òàê êàê äèíàìè÷åñêèé ìàññèâ ýòî óêàçàòåëü
      SyntaxAnalyser.Variable := Variables;
      k1 := fStep * SyntaxAnalyser.calc(err).re;
      if fSolverType = stXY then
      begin
        Variables[0].re := fSolveToResult.Items[i-1].Time + 0.5*step; Variables[0].im := 0;
        Variables[1].re := fSolveToResult.Items[i-1].Value + 0.5*k1; Variables[1].im := 0;
      end else
      begin
        Variables[0].re := fSolveToResult.Items[i-1].Value + 0.5*k1; Variables[0].im := 0;
      end;
      k2 := fStep* SyntaxAnalyser.calc(err).re;
      if fSolverType = stXY then
      begin
        Variables[0].re := fSolveToResult.Items[i-1].Time + 0.5*step; Variables[0].im := 0;
        Variables[1].re := fSolveToResult.Items[i-1].Value + 0.5*k2; Variables[1].im := 0;
      end else
      begin
        Variables[0].re := fSolveToResult.Items[i-1].Value + 0.5*k1; Variables[0].im := 0;
      end;
      k3 := fStep* SyntaxAnalyser.calc(err).re;
      if fSolverType = stXY then
      begin
        Variables[0].re := fSolveToResult.Items[i-1].Time + step; Variables[0].im := 0;
        Variables[1].re := fSolveToResult.Items[i-1].Value + k3; Variables[1].im := 0;
      end else
      begin
        Variables[0].re := fSolveToResult.Items[i-1].Value + k3; Variables[0].im := 0;
      end;
      k4 := fStep* SyntaxAnalyser.calc(err).re;
      Point.Time := i*fStep;
      Point.Value := fSolveToResult[i-1].Value + ((k1+2*k2+2*k3+k4)/6);
      fSolveToResult.Add(Point);
      inc(i);
    until fSolveToResult[i - 1].Value > Limit - 1.2;
  finally
    SyntaxAnalyser.Free;
  end;

end;

{ TProcessRegion }

function TProcessRegion.Add(Item: TDynPoint): Integer;
var
  Temp: PDynPoint;
begin
  New(Temp);
  Temp^ := Item;
  result := inherited Add(Temp);
end;

procedure TProcessRegion.Delete(Index: Integer);
begin
  Dispose(PDynPoint(inherited Get(Index)));
  inherited Delete(Index)
end;

function TProcessRegion.Get(Index: Integer): TDynPoint;
begin
  result := PDynPoint(inherited Get(index))^;
end;

procedure TProcessRegion.Put(Index: Integer; const Value: TDynPoint);
begin
  PDynPoint(inherited Get(index))^ := Value;
end;

function CopyFunction(const Source: TFunction): TFunction;
var
  i: integer;
begin
  SetLength(Result[0], Length(Source[0]));
  SetLength(Result[1], Length(Source[1]));
  for i := 0 to high(result[0]) do
  begin
    result[0, i] := Source[0, i];
    result[1, i] := Source[1, i];
  end;
end;

end.

